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Abstract 

An extended version of the Monte Carlo program GALUGA is presented for 
the computation of two-photon production in e + e~ collisions. Functions im- 
plemented for the five 7*7* structure functions now include several ansatze 
of the total hadronic cross section based on the BFKL-Pomeron and various 
Regge-like models. In addition, structure functions for resonance formation 
are included with full dependence on the two photon virtualities Q\ and Q\ 
as given in the constituent-quark model. The six lowest-lying resonances 
of each of the C-even mesons with J p = 0~, + , 1 + , 2 + and 2~ are pro- 
vided. The program can also be used to calculate with exact kinematics 
the effective two-photon luminosity function. Special emphasis is put on a 
numerically stable evaluation of all variables over the full Q\ range while 
keeping all dependences on the electron mass and Q\. 
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Program Summary 



Title of program: 
Program obtainable from: 



GALUGA 

G.A. Schuler, CERN-TH, 
CH-1211 Geneva 23, Switzerland; 
Ger har d . S chuler @ cern . ch 



Licensing provisions: 

Computer for which the program is designed and 
others on which it is operable: 
Operating system under which the program has 
been tested: 

Programming language used: 
Number of lines: 
Keywords: 



Monte Carlo, two-photon, e + e , az- 
imuthal dependence 
VEGAS (included, 229 lines) 
RANLUX § (included, 305 lines) 
HBO OK § and DATIME J§ for 
the test program (365 lines) 



none 



all computers 



UNIX 



FORTRAN 77 
2209 



Subprograms used: 



Nature of physical problem: 

Hadronic two-photon reactions in a new energy domain are becoming accessible with 
LEP2. Unlike purely electroweak processes, hadronic processes contain dominant non- 
perturbative components parametrized by suitable structure functions, which are func- 
tions of the two-photon invariant mass W and the photon virtualities Q\ and Q2- It is 
hence advantageous to have a Monte Carlo program that can generate events with the 
possibility to keep W and, optionally, Qi at fixed, user-defined values. Moreover, at least 
one program with an exact treatment of both the kinematics and the dynamics over the 
whole range m 2 ^> m 2 (PU/y / i) 4 ^ Qf ~ s (m is the electron mass and yfs the e + e~ cm. 
energy) is needed, (i) to check the various approximations used in other programs, and 
(ii) to be able to explore additional information on the hadronic physics, e.g. coded in 
azimuthal dependences. 
Method of solution: 

The differential cross section for e + e~ — > e + e~X at given two-photon invariant mass W is 
rewritten in terms of four invariants with the photon virtualities Qi as the two outermost 
integration variables (next to W), in order to simultaneously cope with antitagged and 
tagged electron modes. Due care is taken of numerically stable expressions while keeping 
all electron-mass and Qi dependences. Special attention is devoted to the azimuthal 
dependences of the cross section. Cuts on the scattered electrons are to a large extent 
incorporated analytically and suitable mappings introduced to deal with the peaking 
structure of the differential cross section. The event generation yields either weighted 
events or unweighted ones (i.e. equally weighted events with weight 1), the latter based 
on the hit-or-miss technique. Optionally, VEGAS can be invoked to (i) obtain an accurate 
estimate of the integrated cross section and (ii) improve the event generation efficiency 
through additional variable mappings provided by the grid information of VEGAS. The 
program is set up so that additional hadronic (or leptonic) reactions can easily be added. 
Typical running time: 

The integration time depends on the required cross-section accuracy and the applied cuts. 
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For instance, 13 seconds on an IBM RS/6000 yields an accuracy of the VEGAS integra- 
tion of about 0.1% for the antitag mode or of about 0.2% for a typical single-tag mode; 
within the same time the error of the simple Monte Carlo integration is about 0.5% for 
either mode. Event generation with or without VEGAS improvement and for either tag 
mode takes about 4 x 10~ 4 (2 x 10~ 3 ) seconds per event for weighted (unweighted) events. 

Differences with earlier version /j^/: 

i) The ^-integrated total e + e~ cross section a can now be calculated besides da/dW 2 
and d 2 a/dW 2 dQ 2 2 . 

ii) The program can be used to calculate a for a total, Regge-like hadronic cross section as 
well as the effective two-photon luminosity function C defined by a = J dr £(r) cr 77 (rs), 
where <t 77 is the real-photon cross section. In both cases four different ansatze for the Qi 
dependence of the hadronic form factors is provided. 

iii) The total two-photon cross section at large Q 2 calculated in perturbative QCD, in 
terms of the BFKL Pomeron, is incorporated. 

iv) Structure functions for the formation of resonances in two-photon collisions are in- 
cluded for 30 mesons, light or heavy. One can choose between two models: in one, the 
full, non-trivial Q\ , Q 2 correlations as given in the constituent quark model are kept. The 
alternative model is based on a factorized VMD-inspired ansatz. 
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1 Introduction 



Two-photon physics is facing a revival with the advent of LEP2. Measurements of two- 
photon processes in a new domain of 77 cm. energies W are ahead of us |J. Any two- 
photon process is, in general, described J?J by five non-trivial structure functions (two more 
for polarized initial electrons). Purely QED (or electroweak) processes are fully calculable 
within perturbation theory. Several sophisticated Monte Carlo event generators exist 
|8|j9|JT0|] to simulate 4-fermion production in e + e~ collisions. Indeed, the differential cross 
section is not explicitly decomposed as an expansion in the five 7*7* structure functions. 
Rather, the full matrix element for the reaction e + e~ — > e + e~£ + £~ is calculated as a 
whole, partly even including QED radiative corrections. Such a procedure is, however, 
not possible for hadronic two-photon reactions since the hadronic behaviour of the photon 
is of non-perturbative origin. The decomposition into the above-mentioned five structure 
functions (and their specification, of course) is hence mandatory for a full description of 
hadronic reactions. 

Monte Carlo event generators for hadronic two-photon processes can be divided into 
two classes. Programs of the first kind (|TI],|T^.|T^.|T^.[T^|T(J put the emphasis on the QCD 
part but are (so far) restricted to the scattering of two real photons. The two-photon 
sub-processes are then embedded in an approximate way in the overall reaction of e + e~ 
collisions. A recent discussion of the so-called equivalent-photon approximation can be 
found in [17]. 

The other type of programs ]TS|P^[ treat the kinematics of the vertex e + e~ — > e + e~77 
more exactly, but they contain only simple models of the hadronic physics. Moreover, 
the event generation is done in the variables that are tailored for ee — > ee'yy, namely the 
energies and angles (or virtualities) of the photons and the azimuthal angle between the 
two lepton-scattering planes in the laboratory system. Hence, both the hadronic energy 
W and the azimuthal angle <p in the photon c.m.s. (which enters the decomposition of 
the e + e~ cross section into the five hadronic structure functions) are highly non-trivial 
functions of these variables^. 

In the study of hadronic physics one prefers to study events at fixed values of W. Not 
only is W the crucial variable that determines the nature of the hadronic physics (total 
7*7* cross section, resonance formation, etc.), but 77 collisions can be compared with 7p 
and pp ones through studies of events at fixed W PD[ . Next to W, the virtualities Q\ and 
Q 2 of the two photons determine the hadronic physics. At fixed values of W and one of 
the Q's, say Qi, one obtains the cross section of deep-inelastic electron-photon scattering. 
Varying Q 2 , one can investigate the so-called target-mass effects, i.e. the influence of non- 
zero values of Q2 on the extraction of the photon structure function F 2 . Hence it is 
desirable to have an event generator that allows keeping W fixed (or integrating over W) 
and in which Qi and Q2 are the (next) outermost integration variables, so that W and/or 
W and Qi can be held constant. 

The remaining two non-trivial integration variables, which complete the phase space 
of e + e~ — > e + e~X, should be chosen such that three conditions are fulfilled. First, cuts on 
the scattered electrons are usually imposed in experimental analyses. Hence, the efficiency 
and accuracy of the program is improved if these can be treated explicitly rather than 
incorporated by a simple rejection of those events that fall outside the allowed region. 



1 The only program that contains the ^-dependences is TWOGAM |18 . However, the expressions 
taken from [Q are numerically very unstable at small Qf, see the discussion following ([l2]). Moreover, 
itself is not calculated. 
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Second, the peaking structure of the differential cross section should be reproduced as 
well as possible in order to reduce the estimated Monte Carlo error and to improve the 
efficiency of the event generation. And third, it should be possible to achieve a numerically 
stable evaluation of all variables needed for a complete event description. These three 
conditions are met to a large extent by the choice of subsystem squared invariant masses 
Si and S2 as integration variables besides Q\ and Q\. In the laboratory frame, S{ are 
related to the photon energies u>i by si/ 2 — m 2 = 2u 2 /i\^s, where m denotes the electron 
mass. 

In the interest of those readers not interested in calculational details, the paper starts 
with a presentation of a few results in section |^. The differential cross section for the 
reaction e + e~ — > e + e~X is rewritten in terms of the four invariants Q 2 and Sj (i = 1,2) 
in section |3] where also models for the cross section (7(7*7* ~~ > X) are described. The 
integration boundaries with Q\ and Q\ as the two outermost integration variables at fixed 
W are specified in section The derivation of the integration limits is standard |H| but 
tedious. Here the emphasis is put on numerically stable expressions^. To our knowledge, 
numerical stable forms of and <fi are presented here for the first time. All dependences on 
the electron mass and the virtualities of the two photons are kept. The formulas are stable 
over the whole range from Q 2 m - m ~ m {W/^fsY < m 2 up to Q 4 2 max ~ s, i.e. the program 
covers smoothly the antitag and tag regions. An equivalent-photon approximation is 
also implemented (section ||). The complete representation of the four- momenta of the 
produced particles in terms of the integration variables is given in section ^. Section [F] 
describes the incorporation of cuts on the scattered electrons. Details of the Monte Carlo 
program GALUGA are given in section H. 



2 A few results 

In order to check GALUGA, we include the production of lepton pairs, for which sev- 
eral well-established Monte Carlo generators [pjP,[T0| exist. The five structure functions 
for 7*7* — > £ + £~ as quoted in |7[] have been implemented. For the comparison we have 
modified the two-photon part of the four-fermion program DIAG36 (i.e. DIAG36 re- 
stricted to the multiperipheral diagrams) in such a way that it can produce events at fixed 
values of W. The agreement is excellent. Two examples are shown in Fig. [I], the first 
corresponding to a no-tag setup and the second to a single tagging mode. 

Next we study the (integrated) total hadronic cross section. Figures |] and |3|-|5] compare 
different ansatze for the Q 2 behaviour of the various cross sections for transverse and lon- 
gitudinal photons. The results of the two models of generalized-vector-meson-dominance 
type (GVMD (|T7|) and VMDc (pT8|) , dash-dotted and dotted histograms, respectively) 
are hardly distinguishable in the no-tag case, but may deviate by more than 20% in a 
single-tag case. In the contrast, the different Q\ behaviour of a simple p-pole (dashed 
histograms) shows up already in the no-tag mode. Note that this model includes scalar 
photon contributions, but does not possess an 1/Q 2 "continuum" term for transverse pho- 
tons. These differences imply that effects of non-zero Q 2 values must not be neglected for 
a precision measurement of a^iW 2 ). 

During the course of the LEP2 workshop, sophisticated programs to generate the 
full (differential) hadronic final state in two-photon collisions have been developed |22 . 



2 A similar phase-space decomposition with s\ replaced by A = [(s — 2m 2 )(W 2 + Q\ + Q\) — {s\ 
Q\ — m 2 )(s2 + Q\ — m 2 )]/4 is presented in 
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The description of hadronic physics with one (or both) photons off-shell by virtualities 
Q 2 <C W 2 is still premature. Indeed, existing programs are thus far for real photons and 
hence use, in one way or another, the equivalent-photon approximation (EPA) to embed 
the two-photon reactions in the e + e~ environment. It is hence indispensable to check the 
uncertainties associated with the EPA. Hadronic physics is under much better theoretical 
control for deep-inelastic scattering, i.e. the setup of one almost real photon probed by 
the other that is off-shell by an amount Q 2 of the order of W 2 . Corresponding event 
generators exist |22] but also in this case it is desirable to check the equivalent-photon 



treatment of the probed photon. 

An improved EPA has recently been suggested in | I7|| . In essence, the prescription 



consists in neglecting Q 2 w.r.t. W 2 in the kinematics but to keep the full Q 2 dependence in 
the 7*7* structure functions. In addition, non-logarithmic terms proportional to m 2 /Q 2 
in the luminosity functions are kept as well. The study |17| shows that this improved EPA 
works rather well for the integrated e + e~ cross section. In Fig. ^ we show that this EPA 
(solid compared to dash-dotted histograms) works well also for differential distributions, 
with the exception of the polar-angle distribution of the hadronic system at large angles, 
where it can, in fact, fail by more than an order of magnitude! (There, of course, the 
cross section is down by several orders.) 

The EPA describes also rather well the dynamics of the scattered electrons in the 
single-tag mode except in the tails of the distributions (Fig. |3|). The same holds for the 
distributions in the photon virtualities, see Fig. f|. Sizeable differences do, however, show 
up (Fig. |J) in the distributions of the subsystem invariant masses ^/sl. These then lead 
to the wrong shapes for the energy and momentum distributions of the hadronic system 
shown in Fig. |5|. The EPA should, therefore, not be used for single-tag studies. 

Finally we study the prospects of a determination of additional structure functions 
besides F 2 . One such possibility was outlined in |J, namely the study of the azimuthal 
dependence in the 77 c.m.s. between the plane of the scattered (tagged) electron and the 
plane spanned by the beam axis and the outgoing muon or jet. Here we propose to study 
the azimuthal angle between the two electron scattering planes, again in the 77 c.m.s. 
Although such a study requires a double-tag setup, the event rates need not be small, 
since one can fully integrate out the hadronic system but for its invariant mass W. In 
order to demonstrate the sensitivity of such a measurement we show, as a preparatory 
exercise, the distribution for muon-pair production in Fig. ^. Fitting to the functional 
form ^ 

— ^ oc 1 + A\ cos + A 2 cos 2 , (1) 
d0 

we find 

Ai = 0.098 , A 2 = -0.028 . (2) 

Let us emphasize that the selected tagging ranges have in no way been optimized for such 
a study. Nonetheless, given the magnitudes of A4, a measurement appears feasible. 

All but one || event generators for two-photon physics use the azimuthal angle be- 
tween the two scattering planes in the laboratory frame as one of the integration variables. 
In fact, appears as a trivial variable in these programs. None of these up to now pro- 
vides the calculation of 0. An expression for in terms of U, 0, and two other invariants 
is given in |7| (see fl55| ) below) and, in principle, is available in TWOGAM |I8"| . However, 
the factor \Jt\t 2 appears explicitly in the denominator of cos but not in its numerator. 
Hence, at small values of — i,- this factor will be the result of the cancellation of several 
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much larger terms, rendering this expression for cos numerically very unstable. (Recall 
that | ti| mini ~ w?(W/^Js) A <C m 2 , while the numerator contains terms of order s.) In 
contrast, we use the numerically stable expression given in (|56|)]. 
An approximation for in terms of is proposed in ||23|| : 



. 2 Qi Qi (2 s — s l ~ s 2 ) ,,- 
cos a pp rox = cos + sin = . (3) 



Indeed, the correlation between and its approximation is very high in the no-tag case, 
where, however, the dependence on is almost trivial (i.e. flat). Figure |6| exhibits that 
there is still a correlation for a double-tag mode, but formula @ fails to reproduce the 
correct dependence: a fit to (P yields A\ = 0.084 and A 2 = 0.017, quite different from 



3 Notation and cross sections 

Consider the reaction 

e + (p a ) + e~(p b ) -> e + ( Pl ) + X(p x ) + e~(p 2 ) (4) 

proceeding through the two-photon process 

l(qi)+l(q 2 ) ^ X(p x ) . (5) 

The cross section for (|4j) depends on six invariants, which we choose to be the e + e~ cm. 
energy i/i, the 77 cm. (or hadronic) energy W, the photon virtualities Qi, and the 
subsystem invariant masses ^/si•. 

S = (Pa+Pb) 2 , W 2 = Px , 
Sl = (Pi + Pxf = (Pa + <?2) 2 , -Ql =tl=q\= (Pa ~ Plf , 

s 2 = (P2 + Pxf = (Pb + qif , -Ql = h = q 2 = (Pb - P2) 2 ■ (6) 
We find it convenient to introduce also the dependent variables: 

u 2 = si - m 2 - t 2 , v = - (W 2 -ti~ t 2 ) , 

" 1 



Ul = s 2 - m 2 - t x , K = — V /A(W 2 , t lt t 2 ) = — ^v 2 - t x t 2 



\ 4m 2 / 4m 2 

where \(x,y, z) = (x — y — z) 2 — Ayz and m denotes the electron mass. Note that K is 
the photon three-momentum in the 77 c.m.s. In terms of these variables the e + e~ cross 
section at fixed values of y/s and r = W 2 /s is given by: 

da[e+e-^e+ + e- + X] _ a 2 KW , _ ^ 2 ^ 2 ~ 2 , 



dr ~2tt 4 Q 2 Q 2 /3 



dR 3 Z(W 2 , Qf, Q 2 , Sl , s 2 , 0; s, m 2 ) , (8) 



3 This form of </> could, with only minor modifications, be implemented in jq]. 
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where R 3 is the phase space for (|). 

We also give the relation between the cross section at fixed values of r and Q\ and 
the usual form used in deep-inelastic scattering: 



da x s da 



drdt 2 Ql dxdQl 
where x is the Bjorken-x variable defined by 



(9) 



x = 9i = 9i do) 

The hadronic physics is fully encoded in five structure functions. Three of these 
can be expressed through the cross sections a a b for scalar (a, b = S) and transverse 
photons (a, b = T) (asr = 0rs(ft <-> ft)). The other two structure functions ttt and 
tts correspond to transitions with spin-flip for each of the photons with total helicity 
conservation. Introducing 0, the angle between the scattering planes of the colliding e + 
and e~ in the photon c.m.s., these structure functions enter the cross sections as: 

S = 2 p++ 2 p++ a TT + 2 p++ p°° a TS + p?° 2 p++ a ST + pf pf a ss 

+ 2 \pf~ pt~\ r TT cos20- 8 |p+° p+° I tt5 cos0 . (11) 

The density matrices of the virtual photons in the 77-helicity basis are given by 

/ \2 A ? 

'Uo — v) , 4 m 



2pt -low- 1 

{u 2 - v 



K 2 W 2 tj 
2 



K 2 W 2 



V(^+i)ipn = w^- 1 (12) 



with analogous formulas for photon 2. 

A few remarks about the numerical stability of the ^-dependent terms are in order. 
Thus far, these terms are implemented solely in the TWOGAM |TB| event generator, 
using the formulas quoted in 0. Given in and coded in |Tj| are the products X 2 = 



2 \pi~ p^,~ I cos 2 and X\ = 8 |p^° p2~°| cos in terms of invariants. Now, the expressions 
for Xi contain explicit factors of tit 2 (A 2 ) and \Jt\t 2 (X\) in the denominators but not 
in the numerators. Clearly, the evaluation of Xi becomes unstable for small values of 
|tj|. On the other hand, the factors multiplying cos0 and cos 20 in Xi approach perfectly 
stable expressions in the limit m 2 /W 2 — > and ti/W 2 — > 0: 



, 2 ^ . 2m 

\pT H- 2 -( 1 - a; i) + — - 



2 



\ P i°\^ 2 —^^\, (13) 

Xi v 

where Xi = W 2 /si ~ Sk/s (i 7^ k). Hence a numerically stable evaluation of guarantees 
a correct evaluation of the 0-dependent terms. 

The structure functions a a b and r a b for lepton-pair production are often quoted in the 
literature; the formulas of [0] are implemented in the program. Much less is known about 
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the structure functions for hadronic processes. Since we are not aware of a model for r a b 
of the total hadronic cross section, the current version of the program assumes 

ttt = = tts ■ (14) 

The program is set up in such a way that it is straightforward to add a model for r a b- For 
resonance production, r a f, as given in the constituent quark model are implemented. 

The Qi dependence of the cross sections o a b reflects the hadronic physics of the process 
under consideration. For the total hadronic cross section, four Regge-based models are 
provided. They are based upon the assumption 

v ab (W 2 , Qi) = K(Qj) h b (Ql) a^(W 2 ) , (15) 

which is valid for Q 2 -C W 2 ; this is justified in most applications. Note the cross section 
for the scattering of two real photons a^iW 2 ) that enters as a multiplicative factor in 
(0). We take it as |2(J 

a yy (W) =Xs t + Ys- r >. (16) 

The program can be used to calculate a two-photon luminosity function if one takes 
a 17 (W) = 1. 

The four models are defined as follows. The first one is based upon a parametrization 
fl24j of the 7*p cross section calculated in a model of generalized vector-meson dominance 
(GVMD): 



h T (Q 2 )=rP 1 - 2 (Q 2 ) + (l-r)P 2 -\Q 2 ) 

h s (Q 2 )=dr%Pr 2 (Q 2 ) + (l-r] 
mf 



lnP 2 (Q 2 )-P 2 -\Q 2 



^ 1 , Q 2 



Q 2 

P t (Q 2 ) = l + ^, (17) 
mf 

where we take £ = 1/4, r = 3/4, mf = 0.54 GeV 2 and m\ = 1.8 GeV 2 . 

The second model |25] adds a continuum contribution to simple (diagonal, three- 
mesons only) vector-meson dominance (VMDc): 



)2 



12 / ^,2 \ 2 



w<n -^§^feft?J ' (18) 

where r p = 0.65, = 0.08, = 0.05, and r c = 1 — XV r v- 

Since photon-virtuality effects are often estimated by using a simple p-pole only, we 
include also the model defined by (p-pole): 

MQ 2 )=(-^V , h s (Q 2 ) = ^(^±-) 2 . (19) 
\m 2 p + Q 2 J m 2 p \m z p + Q z J 

The fourth model is identical to flT§|) but has hs(Q 2 ) = 0. 
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i\ j 




1 


2 


3 


4 


5 




n 2S+1 Lj 


J = 1 


1 = 


(/ = 0)' 


cc 


bb 


1 




7T° 


n 


rf 




\m (9400)1 

Li" \ /J 


2 


1 3 P 


a (980)* 


f (980)* 


/o(1370)* 




A- OO 


3 


1 3 P 1 


ai(1260) 

J- V / 


A(1285) 


/i(1510) 




/vox 


4 


1 3 P 2 


a 2 (1320) 


/ 2 (1270) 


/'(1525) 






5 


1^2 


vr 2 (1670) 


[^(1680)] 


[^(1890)] 


lW3840)] 


K (10150)] 


6 


2% 


tt(1300) 


77(1295) 


[77 r (1400)] 


r/ c (2S) 


M9980)] 



Table 1: Lowest-lying C = +1 mesons an d their labelling with i and j. States in 
brackets are not yet found. Status of states marked with a * is not yet clarified. 



At large virtualities the behaviour of the 7*7* cross sections is fully predicted by 



perturbative QCD in terms of the BFKL Pomeron | 2q| . We use the results obtained in 
the so-called saddle-point approximation 



<Tab = W a W { 

where W T = 9, W L 

N c a 



7T 




exp (4£ In 2) 



exp 



56^(3) 



(20) 



a. 



(33-2n/) ln(/i 2 /A 



Q = c q Qi Qi , c q 
2 - CuQi Q2 



10 2 



exp(-5/3)) 



(21) 



In order to ensure the validity of the high-energy approximation that went into the cal- 
culation of |2f| we demand 



Qm'm ^ Qi ^ Qr 



W 2 >SQ 1 Q 2 (5=10 2 



(22) 



Structure functions for resonance formation in two-photon fusion were recently cal- 
culated in the constituent- quark model [p7| . Although the results stricly apply to heavy 
mesons only, the Qi dependence is presumably also very reasonable for the lighter mesons. 
The mesons included are listed in Table [I]. The structure functions are given by 



a ah [r\ = (2J + 1)8tt 2 ^ f ab (Q 1 ,Q 2 ,M p ; J F )BW(W,M) . 



(23) 



Here M denotes the mass, J the total spin, P the parity, and T the total width of the 
C = +1 meson. The mass M p is equal to M for all mesons except for it , 77, and rf, for 
which we take the p mass. T 77 is the two-photon decay width for all mesons except for 
those with J = 1, where a different quantity had to be introduced since J = 1 mesons 
cannot decay into two real photons. Explicit expressions for T 77 and fab can be found in 
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|2T|| . Form factors for the interference terms t?t and Txs are also implemented. The W 

BW (W, M) 



dependence is given by 



7T (m 2 - w 2 ) 2 + M 2 r 2 

= 5(W 2 - M 2 ) , (24) 

depending on whether one integrates over W or keeps W fixed. 

Note that the form factors f a b(Qi, Q2, M p ; J p ) do not factor in Qi- and (^-dependent 
factors, nor do they have simple monopole or dipole behaviours. As an alternative a 
simple factorizing model based on VMD is also implemented 

MJ P \ - & + (^) 2 (A) 2 BWffiM) . (25) 

Observe that all 1 + cross sections are zero for model (ET~ 



4 Phase space 

The phase space can be expressed in terms of four invariants^]: 

= ttit I dt2 dtl dsi dS2 -Af= ' ( 26 ) 

lop s J v — A 4 

where A 4 is the 4x4 symmetric Gram determinant of any four independent vectors 
formed out of p a , p&, pi, px, Pi- The physical region in t 2 , ti, si, s 2 for fixed s satisfies 
A4 < 0. Since A4 is a quadratic polynomial in any of its arguments, the boundary of the 
physical region, A 4 = 0, is a quadratic equation and has two solutions. Picking s 2 as the 
innermost integration variable, the explicit evaluation of A4 yields 

16 A 4 = a s\ + b s 2 + c = a (s 2 - s 2+ ) (s 2 - s 2 _) , (27) 

where 

a = X(si,t 2j m 2 ) 

b=-2s m 2 t x - 2 m 2 s\ + 8t 2 m 4 - 2 m 2 t\ - 2 s si W 2 + 2 m 2 s W 2 + 2 ij s s 1 + 2 s t 2 a x 
+ 4 m 2 si W 2 + 4 m 4 si + 2 t x t 2 s - 2 i 2 m 2 ii - 2 t 2 s - 2 m 2 t 2 s + 2 £1 1 2 «i 

- 4 m 4 iy 2 - 2 t x s\ + 2 s t 2 W 2 - 2 m 6 + 2 m 4 ti 

c = -2 s m 4 W 2 - 2 t 2 m 2 si - 2 t x i 2 s 2 + 2 s £1 t 2 s x - 2 s t\si + £ 2 s 2 + t 2 s 2 + t\s 2 

+ m 4 s 2 + mH\ - 6 m 6 ti - 2 m 6 si - 4 m 4 si W 2 + 2 m 4 £ 2 s + 2 m 4 i 2 s x + 8 m% si 

- 2 s 2 £ 2 W 7 2 - 2 tx s 2 M/ 2 - 2 m\ s\ + m 8 - 2 m 2 s £ 2 «i + 4 m 6 jy 2 + m 4 £ 2 
+ 4 m 2 ti t 2 s - 2 m 2 ti £ 2 «i - 6 m% + s 2 W 4 + 6 m 2 s t 2 W 2 -As m 2 W 4 
-2s m 2 t\ + 2 s ti m 4 - 2 s m^ 2 . + 2 ti t 2 m 4 + 2 s m 2 Si W 2 - 4 s h t 2 W 2 

-2 s t Y m 2 Sl + 6st t m 2 W 2 + 2st lSl W 2 . (28) 



4 For the fully differential cross section a factor 2ir has to be replaced by a trivial azimuthal integration 
around the z-axis. 
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A numerical stable form for the s 2 limits is 

-b+ y/E 
S 2+ = o 

2 a 

S2- = — , (29) 
as 2+ 

where A = b 2 — Aac is given below in a numerically stable form, in (|32|). 

In order to remove the singularity due to (— A4) -1 / 2 (in the limit m 2 <C Sj, VT 2 , 
the s 2 integration degenerates to an integration over the <5-function S(s 2 — sW 2 /si)), it 
is advisable to change variable from s 2 to x&, < x& < 1: 

s 2 = — {-6 - cos(a;4 7r)| 

ds 9 4 7r r 1 . , 

-= / dx 4 . (30) 



82- V~^4 

For later use we also need a numerically stable form of the Gram determinant, which 
reads 

16A4 = - Asi " 2(X47r) . (31) 
4a v ' 

The Si-integration limits follow from the requirement A > 0. They are most easily 

derived when realizing that the discriminant A is given as the product of two 3x3 

symmetric Gram determinants or, equivalently, the product of two kinematic G functions 



^A = 4G 3 G 4 = 64L> 3 L> 4 > ( 32 ) 



where 



- 4 D 3 = -4 A 3 (p a , p b , q 2 ) = G(s, t 2 , St, m 2 , m 2 , m 2 ) = G 3 

-4 D 4 = -4 A 3 (p a , q u q 2 ) = G(t h 81, t 2 , m 2 , m 2 , W 2 ) = G A . (33) 

Since any 3x3 Gram determinant A 3 satisfies A 3 > 0, the physical region is that where 
both G 3 and G4 are simultaneously negative. Solving Gi for s\ 

G 3 = m 2 (si - s u +) (si - su_) 

= m 2 s\ — 2 m 4 S! -st 2 Si-3 m 2 t 2 s + m 6 + t 2 s 2 + t^s 
G 4 = £1 (si - S12+) (si - S12-) 

= -2 ti m 2 si - t 2 m 2 *i + ^1 - m 2 W% + m 2 t\ + t 2 W 2 ti - t x s x W 2 



v 2 

2 m 2 t 2 W 2 + m 2 Vr 4 + t x s\ + t 2 ^ - t x t 2 «i (34) 



we find 



Sllzt 



t 2 S + 2 m A ± V /A(5, m 2 , m 2 )A(t 2 , m 2 , m 2 ) 
2 m 2 



t 2 2 l^ 2 tx JXitut^W^Xitum 2 ,™ 2 
si2± = — + m + — 7T =•= 



2 2 2 2*! 

t 2 5(-3m 2 + 5 + t 2 ) 



Sn+Sii- = = hm 



■i 



m 2 



S12+S12- 



[W 2 - m 2 ) (-m 2 + t 2 ) + ^ ^- . (35) 
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Note that Si 2+ < s 12 -. Since G 3 is always negative between its two roots, the range of 
integration over s\ is S12- < s\ < sn+. Numerically it is more advantageous to calculate 
the limits as 

simin = s 12 - = m 2 + - (w 2 -t 1 +t 2 + y x ^\{W 2 MM) 
2 2 (s + t 2 - 4 m 2 ) 

Slmax = Sn+ = m H — — . (36) 

1 + P J/a 

The dominant behaviour of the si integration is given by the factor A -1 ' 2 (si, £2; J^ 2 )) 
see ( |30|) . (In the limit t 2 , m 2 < si, this becomes dsi/si integration.) This factor can be 
transformed away by the variable transformation from s x to X3, < X3 < 1, 

si = Xi/2 + m 2 + t 2 + 2 m 2 t 2 /X 1 
X 1 = (u + KW)(l + y 1 ) e W (S 1 x 3 ) 
A -1 s(l + /3) 2 

1_ (- + ^^)(l + 2/l )(l + l/ 2 ) ' ( } 

such that 

dsa — = 4vt(5i / dx 3 . (38) 

Slmin JO 



The physical region in the ti-t 2 plane is defined by the requirement G{ < for all s± 
values between the limits (m + W) 2 < Si < (-y/s — m) 2 . Since for the reaction considered 
here the masses of the particles involved are such that the values t\ = (m a — mi) 2 , 
t 2 = (m& — m 2 ) 2 cannot be reached and t 2 is never larger than zero, the boundary curve 
in the t x -t 2 plane is simply given by si 2 _ = %+. Equivalently, the t x limits can be found 
by solving G4 = with s\ = sn+ for t\\ 

^lrnin = — 7: ( hAtij , tlmax = 7 5 (39) 

* \ a l J a l Hmin 



where 



ai 

r2\ 



a x = 2 (Q + t 2 + 2 m 2 + W 2 
b x = Q 2 - W A + 2 F/ 2 t 2 - t\ - 8 m 2 t 2 - 8 m 2 jy 2 
Cl = 4 m 2 (W 2 - t 2 ) 2 
Ai = b\ - 4aici 

= (Q + t 2 - P7 2 + 4 m W) (Q + t 2 - W 2 - 4 m W) 

(Q 2 - 2 Q t 2 + 2 Q W 2 + 1\ + W 4 - 16 m 2 t 2 - 2 W 2 t 2 ) 

Q = — - \t 2 s — m 2 t 2 — m 2 W 2 + J \(s, m 2 , m 2 )X(t 2 , m 2 , m 2 )\ 
mr I v J 

= 4(. + *,-4m') ^ 

1 + ^2/2 

Finally, the ^-integration limits follow from requiring Ai > 0: 

(^2 — ^21+) (^2 — ^21-) (^2 — ^22+) (#2 — ^22-) t 2 (t 2 — t 23 ) 2 

Al ~~ P P P ' ^ 

Pi t 2 ^3 
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where 

F 2 



4s 

h sfi y 2 + t 2 s - 2 W 2 m 2 + 4 m 3 W 
4s 



h sf3 y 2 + t 2 s - 2 W 2 m 2 - 4 m 3 W 
F 3 = -16 m 4 / { -2 i 2 s 2 /3 y 2 + 4 t 2 s/3 y 2 m 2 -t 2 s 2 - t 2 s 2 (3 2 + 4 1 2 sm 2 

+4 s 2 l3 2 m 2 - 4 t 2 m 4 + 16 m 6 } 



2 -W 2 + s-Am 2 ±f3 Ms - W 2 ) (s - W 2 ) 



<21± 



^22± — 



-2 mPy - W 2 + s - 4 m 2 ± (5 Ms - W 2 ) (s - Wl) 



2^2 



2.3 



(s-2m 2 ) 
m 2 



W± = W ±2m . (42) 
Equivalently, they are arrived at by solving G 3 = with s± — (m + W) 2 for t 2 : 

t 2 min = t 22+ = - X - (s - W 2 - 2 m W - 4 m 2 + Af 2 ) 
m 2 W 2 W 2 

^2max — t 22 - — T , (43) 

S t-2 min 

where 



At 2 = /ty(s-W 2 )(s-W^) . (44) 
The phase space finally becomes 

di? 3 = — 5- / dt 2 / d*i 5i(*i,t 2 ) / drr 3 / dx 4 . (45) 

4pS Jt2min "^*lmin •'0 ^0 

The dominant ^ behaviour is taken into account through a logarithmic mapping, so that 
we end up with a cross section of the form 

% = n / d ^ *w 

TT /dx, ln^ ln^ In- S ?, + ®\, r E ■ (46) 



i=l 



Finally, the total cros section is obtained by integrating over W. The kinematical 
limits are m n < W < Ms — 2m. In the case of resonance formation, a Breit-Wigner 
mapping is performed, while a logarithmic mapping is used for all other cases. 



5 Equivalent-photon approximation 

An approximation is arrived at by neglecting as much as possible the electron-mass and 
ti dependences in the kinematics, but keeping the full dependence on W and Qi in the 
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hadronic cross sections cr ab (W 2 , Q\, Q|) |I7) : 



da _ /■« dsi /■*!» dt 2 dti y* / sW/ 2 \ a 2 W/ 2 

dr Ay 2 s\ Jt 2 a t 2 Jt la h Jw 2 2 \ sj / 16 7r 2 s 

{ 2 Plapprox 2 Plapprox a TT + 2 Pi" a p prax Plapprox a TS 
"^Plapprox 2 Plapprox °"ST + Plapprox Plapprox °~Ss] • (47) 

The integration limits are given by: 

Tfl X; , \ • 2 ^imax 

(1 — Xi) sin 



\ — Xi 2 

t* = - (1 - a*) sin 2 ^ , (46 
1 — x,- 2 



where xi = S2/S and £2 = si/s. 

The approximate forms of the photon density matrices read: 

2 Plapprox ="2 l + (l-^l) 2 - 1 



approx ^2 I 1 \ L > Q2 

4 



Plapprox 2 (1 "^l) ' (49) 



X 



6 Momenta 

Here we present the particle momenta in the laboratory frame. The particle energies 
follow simply from E^ = (p a + p b ) ■ pij ^Js: 



s + m 2 — s 2 



2^S~ 



s + m 2 - Si 
E 2 = 1= 

si + so — 2 m 2 . . 

g * = 2V i (50) 

and the moduli of the three- momenta from P 2 = E 2 — m 2 . The polar angles 0j with 
respect to the beam axis could be calculated from p b ■ Pi = E b E± — P b Pi cos 6i 

s - s 2 + 2t x - 3m 2 
costal = - 



cos 8 



2(3^s~P x 
s — Si + 2t 2 — 3 m 2 



2 



2 



s 2 -s 1 + 2(t 2 -t 1 ) 
C ° S6X = 2j^P~x • (51) 

Typically, the polar angles are very small and it is better to calculate them in a numerically 
stable form from 

■ a 2 ^ 
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" s(3P x ' 

Equations ([y]) are then only used to resolve the ambiguity $i <-» n — Oi. The quantity D 3 
is defined in (|33|-p5D; D 1 is obtained from D 3 by the interchange t 1 <-> t 2 and Sj <-> s 2 . 



The same interchange relates Z?2 ; needed below, with D 4 , given in (|33H35|). Furthermore, 
we have: 



D 5 = D 1 + D 3 



2D 6 

,2\ /tt/2 



7» 



(a - 4m 2 ) (W^ 2 - t x - t 2 ) + (si -t 2 - m 2 ) (s 2 - h 

2 



m 



4 



si 



m 2 ) (s 2 — m 2 ) 



(53) 



The polar angles (f>i (<p 2 ) between the e + (e~) plane and the hadronic plane and the 
polar angle 4> = 4>i + 4>2 between the two lepton planes in the e + e~ c.m.s. are again best 
calculated using the numerically more stable form for the sinus function 



cos< 



sin . 



sin <px 



sin <p 2 



cos (pi 



s/3 v^a; 

s (3 Px sin 9x Pi sin 0\ 

s j3 Px sin Ox P2 sin 2 
Dx + D 6 



COS 02 



r D 3 + \[D\ cos * 



=• (54) 

'D 3 + D t + 2 y^DTD^ cos 

An expression for the azimuthal angle between the lepton planes in the 77 c.m.s. can 
be deduced from the formulas given in @ : 

-2 s + ui + u 2 - u + Am 2 + v(u 2 - v) (u t - v)j{K 2 W 2 ) 



cos 1 



kxt 2 (2p++-2) (2p++-2 
Numerically more stable is the following form 



(55) 



smi 



KW V=A 4 



and 



COS 1 



where 

\<6D 7 = 2 W 2 (s x s 2 



sW 2 - 2t x [-t 1 s 1 + st x + sx s s 



W 2 Sl -2sW 2 



(56) 



2t 2 i-t 2 s 2 + t 2 s + Sl s 2 -2sW 2 + s 2 W 2 



2txt 2 [sx + 2s + 2W 



S-2 



-2m 2 (m 2 t 2 - t 2 2 - m 2 W 2 + m 2 t x -2W i -t 1 2 + W 2 s 1 + 2 M 2 



+ 3ti W 2 -txSx + 3t 2 W' 



t 2 s 2 — to. s 



2 »1 



s 2 W 2 -t x s 2 



(57) 
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A numerically stable relation between and at — £,, m? <C W 2 is provided by 

4 s 2 = 4 $2 -D3 + 2 £2 s 2 r2 cos <p s 2 — 2 £x £2 SS2 (s — Si) 

— sti £ 2 ( _ 2 £ x £ 2 — sti + ti Si + 3 £ 2 s 2 — t 2 s + 2 r 2 cos s) 

— Am sr 2 cos0 Si s 2 + O (m s x s 2 /s, m U s Si s 2 J , (5? 

an analogous expression for .D 4 , and 



16 s J D 2 D± cos = 16 W 2 JD 1 D 3 cos — 4 1 2 s 2 r 2 cos — 4 £x s 2 r 2 cos 

+ 4 ti £ 2 s (— si + 2 s - s 2 + U + t 2 ) + 8 m 2 cos r 2 si s 2 
2m 2 £ 2 



Z /At to / 9 9 

H -t 2 s + 4 s£ 2 s 2 - 2 £ 2 s 2 - 4 ssi s 2 + 2 si 



S2 



1 r 2 cos ss 2 + s 2 Si + s 2 s 2 + 8 r 2 cos s 2 



+ 



2m 2 £i 



-£x s 2 + 4 £x ssx — 2 ti Si 2 — 4 ssx s 2 + 2 Si 2 s 2 



r 2 cos (f) ss\ + s 2 Si + s 2 s 2 



+ 8 r 2 cos s • 
+ O (m 4 s 2 Sj/s, m 2 ti t 2 s^j 



where 



m 



+ h 1 + -- 



m 



S2 



+ ti 1 



fa 
s 



For m — > and ti/W 2 — > 0, (|58| ) and ( [59"D lead to the approximate relation 
The four-momenta are now given by 

p a = ^(1,0,0,-/3) 
p b =iv^ (1,0,0,/?) 

Px = (-Ei, —-Pi sin^x cos0x, — P± sin^x sin0x, — Px cos#x) 
p 2 = (E 2 , -P 2 sin 9 2 cos0 2 ,P2 sin ^2 sin0 2 ,P2 cos# 2 ) 
p x = (E x ,Px sin6 x ,0,Px cos9 x ) . 

Finally, a random azimuthal rotation around the z-axis is performed. 



(59) 



(60) 



(61) 



7 Experimental cuts 

If cuts on the angle 9 2 and the energy E 2 of the scattered electron are applied, the (si, £2) 
integration region shrinks as follows (see Fig. |7]): 



and 



sxiow = min < (m + W 7 ) , m + s 1 



■Siupp = max < (-y/s - m) 2 , m 2 + s 1 



^M-SlAmax) < £2 < T 2 (si,6> 2 

16 



2E 



2 max 



2 -E2 m in 



(62) 
(63) 



where 

T 2 (s 1 , 9 2 ) = ~ [3 m 2 - s + s 1 + (3 cos0 2 \A(s, si, m 2 ) 



2m 2 ( Sl -m 2 ) 2 1/2 2 . 2 #2 

p A 7 (s, si, m ) sin 



s [/3A 1 /2( S)Sl)m 2) + s _ Si _ 3m 2] ^ v , 1, ; 2 

._{^ +s(1 _ l2) w|}. (64) 

The approximate form holds for m 2 <C Si and a small angle #2 and is used in (^8[) . 
If, as in our case, £ 2 is the outer integration, then its lower limit becomes 

hmin = min {T 2 ( 

Slupp j "2 max Sllow; "2 max )} , (65) 

while the upper limit is more complicated 

^2max T 2 (-Sluppj ^2 min) > Si U pp 

= T 2 ( 

Sllowj ^2 min) Si 
= *2 Snow < Si < Si U pp , (66) 



where 



2 2m y/s 
s\ = s + m — 



m 2 + 



X 

sf3 2 sin 2 9 2 



X (l+2m/VsI 
i 2 = 2 m 2 — m ViX (0 2 < tt/2) 

= 2m — mys t= {" 2 >^/2) 

vX 



I = — + /? 2 siii 2 e 2 . (67) 

s 

The Si-integration range is a rather complicated function of t 2 and may even consist of 
two separated ranges (Fig. |7|). Moreover, the Si-integration range is affected by t\ and 
cuts on Ei and 6\. Then it is better to use the Monte Carlo method. In any case, since 
the ti integration are the most singular ones, the most important constraints are taken 
into account through ( |65D and (|66|) and the analogous formulas for t\. 



8 Details of the program 
8.1 Common blocks 

The user can decide whether to calculate (i) the fully integrated cross section a(s) (in 
Wmin < W < W max and t2min < t 2 < t 2max ) , (ii) the cross section at fixed r, dcr(r; s)/dr 
(|) with W as given in the ggLcrs call (W = M R for resonance production) and t2min < 
t 2 < ^2max, or (iii) d 2 o"(r, £ 2 ; s/dr dt 2 (^) with W as given in the ggLcrs call (W = M R for 
resonance production) and t 2 at the user-defined value t2user. In the case of da/dr, the 
user can choose between the exact or an approximate treatment ( pTTD of the kinematics. 
If lower and upper integration limits lie outside the physical range (H^min > Tn n , W max < 
y/s — 2m, t 2m[n > —s, t 2max £ 0), the full phase space is taken. 
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Common /ggLapp/Wmin,Wmax, t2user ,t2umin, t2umax, iapprx, ivegas , iwaght 



Wmin 

Wmax 

t2user 

t2umin 

t2umax 

iapprx 



Minimum hadronic mass W for iapprx = -1 (Default (D): m % ) 



ivegas 
iwaght 



Maximum hadronic mass W for iapprx = -1 (D: */s — 2m). 
Fixed value of t% chosen by user for iapprx = 1 (D: — 5GeV 2 ). 
Minimum value of t 2 for iapprx ^ 1 (D: — s). 
Maximum value of t 2 for iapprx 7^ 1 (D: 0). 
= — 1: Total cross section integrated over 

Wmin < W < iy max and t 2 min < h < t 2max ; 

d 2 cr/drdi 2 at W as specified in ggLcrs or W = Mr and t 2 = t2user; 
dcr/dr at W as specified in ggLcrs or W = Mr for t2min < h < fomax] 
as iapprx=0 but using approximate kinematics. 
VEGAS integration; 
Simple integration. 
Unweighted events, i.e. Weight = 1; 
Weighted events. 



Cuts on the scattered leptons are set in 

Common /ggLcut/thlmin , thlmax , Elmin , Elmax , th2min , th2max , E2min , E2max 
thlmin,thlmax 



th2min,th2max 



Elmin, Elmax 
E2min,E2max 



Minimum and maximum scattering angles of scattered e" 1 
w.r.t. direction of incident e + . 

Minimum and maximum scattering angles of scattered e" 

w.r.t. direction of incident e~. 

Tighter cuts should be applied to the e~. 

Minimum and maximum energies of scattered e + . 

Minimum and maximum energies of scattered e - . 



Models for the 7*7* cross sections and their parameters are chosen in 



Common /ggLmod/ imodel 



imodel 
imodel 
imodel 
imodel 
imodel 
imodel 
imodel 
imodel 
imodel 
imodel 
imodel 
imodel 



1 

2 
3 
4 
9 

30 
31 
32 
33 
34 

100 + 10z+j 
200 + lOi + j 



for luminosity function (cr 



77 



0V 



= 1); 
i); 



GVMD model 

VMDc model ([18]) for luminosity function v — 
p-pole model flTP| ) for luminosity function (cr 77 : 
as 3, with h s (Q 2 ) = 0; 
Exact cross section for lepton-pair production; 
BFKL model (0) of a*£; 
GVMD model QT7J) for a l a f with a 77 of (]T§); 
VMDc model (0) for a™ with a 77 of (|T§); 
p-pole model fll9D for 0*$ with (J 77 of (|16|); 
as 33, with h s (Q 2 ) = 0; 

Meson cross section (|23|) with i,j according to Table ^ 
Meson cross section (|25| ) with i,j according to Table |I[ 



Common /ggLhad/ r,xi,mls 5 m2s,rrho,romeg,rphi,rc,mrhos. 
& momegs ,mphis ,mzeros 

Parameters for (|I7|-4T9"|) : r, ^, mf, m 2 ,, r p , r w , r^, r c , m 2 p , rri^, rri^, m§. 



18 



Common /ggLres/ Rmass ,Rwidth,Pmass ,Rtotw, UP, iq, il 

Parameters for (^3|) : M, T 77 , M p , T, i, j, int (imodel/100). 

Common /ggBFKL/ delta, Qmin, Qmax, Lambda, Nf , cQ, emu 

Parameters for (gty: 5, Qmin, Qmax, A, n/, cq, c m . 

The integration variables and the particle momenta are stored in 

Common /ggLvar/ yar(10) , 
& t2,tl,sl,s2,El,E2,EX,Pl,P2,PX,thl,th2,thX,phil,phi2,phi,pht 

yar(i) Integration variables for VEGAS. 

t2,tl,sl,s2 Invariants t 2 , t\, Si, S2. 

E 1 , E2 , EX Energies E h E 2 , E X - 

P 1 , P2 , PX Three-momenta P x , P 2 , P x . 

thl,th2,thX Polar angles #1, 9 2 , 6 X . 

phil ,phi2,phi,pht Azimuthal angles (pi, <p2, <P, <P- 

Common /ggLvec/ mntum(7,5) 

Particle four-momenta mntum(i,k): k = 1...5 for p x , p y , p z , E, sign(jo 2 ) x J\p 2 \ 
i = 1. . .7 for incident e + , incident e~, photon from e + , photon from e~, scattered e + 
scattered e~, hadronic system X. 



Parameter for the simple integration and results of the integration and event generation 
are stored in 

Common /ggLuno/ cross, error, Fmax,Fmin, Weight, npts ,nzero ,ntrial 



cross 
error 
Fmax 

Fmin 
Weight 
npts 
nzero 



ntrail 



Estimate of luminosity. 

Estimate of error on luminosity. 

Maximum function value, calculated in ggLcrs; 

checked in ggLgen. 

Minimum function value, calculated in ggLuF . 
Weight if weighted events requested. 

Number of function evaluations for simple integration (D: 10 6 ). 
Number of cases where function was put to zero in ggLuF 
because it failed the cuts; 
initialized to zero in ggLcrs, ggLgen. 

Number of trials necessary in ggLgen to generate an event; 
incremented by each call. 



Parameters for the VEGAS integration are set in 

Common /ggLvgl/ xl(10) ,xu(10) ,acc,ndim,nf call,itmx,nprn 
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acc VEGAS accuracy (Default (D): 1(T 4 ). 

ndim Number of integration variables (D: 4). 

nf call Maximum number of function calls per iteration for VEGAS (D: 10 5 ). 

itmx Number of iterations for VEGAS (D: 4). 

nprn Print flag for VEGAS (D: 2). 

Additional common blocks 

Common /ggLprm/ s,roots,Whad,m,Pi,alem 

s Overall cm. energy square s (D: 10 4 GeV 2 ). 

roots Overall cm. energy ^/s (twice the beam energy), 

set by user through call to ggLcrs . 
Whad Hadronic mass W, set by user through call to ggLcrs (D: lOGeV). 
m Electron mass (D: 511 keV). 

Pi 7T 

alem a em (D: 1/137). 

Common/ggLvg2/XI (50 , 10) , SI , SI2 , SWGT , SCHI , NDO , IT 
Common /ggLerr/ 

& itl , iDl , iD3, iD5 , itX, iph, ipl , ip2, ial , ia2 , ia3, ia4, iel , ie2 , ipt , is 
Block Data ggLblk 

8.2 Subroutines 

ggLcrs (rs,W) Calculates a, da/dr, or dcr/drd^ (see iapprx) 
and finds Fmax; rs = y^s, W = W. 

ggLmom Builds up four-momenta. 

ggLprt Prints four-momenta and checks momentum sum. 

ggLgen(Flag) Generates one event; 

Flag=F if a new maximum is found; then it is advisable 
to restart event generation with adjusted maximum. 

InitMassWidthQ , j ,M,G,GT,PM,pn) Initializes resonance parameters. 

8.3 Double-precision functions 

ggLint(W2,m2,Qls,Q2s,sl,s2,phi,s) E as defined in (|TT|). 



20 



ggLuF(xar ,wgt) 

ggLriTT(W2,Qls,Q2s) 

ggLhTS(W2,Qls,Q2s) 

ggLhSS(W2,Qls,Q2s) 

ggLrTS(W2,Qls,Q2s) 

ggLrTT(W2,Qls,Q2s) 

ggLhT(Qs) 

ggLhS(Qs) 

ggLgg(W2) 

ggLuG(z) 

mucrss (tl ,t2 , i) 

SBFKL(Ql,Q2,i) 

resTT(W2,Qls,Q2s,i) 

resTS(W2,qis,Q2s,i) 

resSS(W2,Qls,Q2s,i) 

tauTT(W2,Qls,Q2s,i) 

tauTS(W2,Qls,Q2s,i) 



defined in (|46|) . 
a TT (W 2 ,Q 2 ,Q 2 ) 
°Ts{W\QlQl) 
°ss(W 2 ,QlQl) 
tts(W 2 ,QIQ 2 2 ) 
ttt(W 2 ,QIQD 
h T (Q 2 ) 
h s (Q 2 ) 

a^iW 2 ) or 1 depending on imodel 

Makes the variable transformation from Xi in ( 46|) 
those used by the simple or VEGAS integration. 

Muon-pair cross sections 

BFKL cross sections 

ott for resonances 

(Tts for resonances 

ass f° r resonances 

ttt for resonances 

tts for resonances 
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8.4 Excerpt from the demonstration program 

* Initialize the random number generator RanLux 

Call rLuxGo (3,314159265,0,0) 

* 

* Initialize GALUGA; get luminosity within cuts 

Call ggLcrs(rs,W) 

* 

* Initialize plotting 

Call User(0) 

* 

* Timing: 

Call Timex(timel) 

Call rLuxGo (3, 314159265, 0,0) 

* 

* Event loop 

Do 10 i=l,Nev 

Call ggLgen(Flag) 

If ( .not .Flag) Write(6,*) 'Caution: new maximum' 

* 

* Calculate 4-momenta 

Call ggLmom 

* 

* Display first 3 events 

If(i.le.3) call ggLprt 

* 

* Fill histrograms 

Call User(l) 
10 Continue 

* 

Call Timex(Time2) 

Write(6,300) Nev,Time2-Timel , (Time2-Timel)/real(Nev) , 
& iwaght,ntrial,nzero,Fmax 

* 

* Finalize plotting 

Call User(-Nev) 

* 

300 Format (/,3x, 'time to generate ',18,' events is ',E12.5,/, 
&3x, 'resulting in an average time per event of ',E12.5,/, 
&3x, 'unweighted events requested if 1 : ',18,/, 
&3x,'the number of trials was: ',18,/, 
&3x,'the number of zero f was: ',18,/, 
&3x,'the (new) maximum f value was: ',E12.5) 

* 

Stop 
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Figure 1: Comparison of muon-pair production in GALUGA (dashed histograms) and 
DIAG36 (solid histograms) at ^/s = 130 GeV and W = 10 GeV. Top: distribution in the 
logarithm of the polar angle of the fi + [i~ system; no cuts are applied on the scattered 
electrons. Bottom: distribution in t\ under the cuts: 1.55 < Q\ < 3.67° and 30 GeV< E\. 
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Figure 2: Distributions in Ei, ln#i, Ex, and 6x for the integrated total hadronic cross 
section at y/s = 130 GeV and W = 10 GeV. No cuts on the scattered electrons are 
applied. Histogram line-styles correspond to GVMD model in the EPA (solid), p-pole 
model (dashed), GVMD model (dash-dotted), VMDc model (dotted). 
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Figure 3: Distributions in E 1: E 2 , ln^, and ln# 2 for the integrated total hadronic cross 
section at y/s = 130 GeV and W = 10 GeV. The cuts 9 1 < 1.43°, 1.55 < 9 2 < 3.67°, and 
30 GeV< E 2 have been applied. Histogram line-styles correspond to GVMD model in the 
EPA (solid), p-pole model (dashed), GVMD model (dash-dotted), VMDc model (dotted). 
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Figure 4: Same as Fig. [3], but for the distributions in y/sl, ln(— s/ti), and £2 
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Figure 6: At the top, the correlation between ap prox ©, proposed in [23|| , and 0; at the 
bottom, the distribution in (solid histogram) and its approximation (dashed histogram) 
for the integrated muon-pair cross section at y/s = 130 GeV and W = 10 GeV. The cuts 
1.55° < 6i, 5 GeV< E±, and 30 GeV< Ei have been applied. Also shown is a fit to dcx/d0 
of the form 1 + A\ cos <p + A2 cos 20. 
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Figure 7: Phase space in the variables (t 2 ,si) for yfs = 4 and m — 1. The solid lines 
correspond to #2 = 0, 7r/4, 7r/2, 37r/4, and 7r (from t 2 = to t 2 = —12 at Si = 1). The 
dashed lines are si = §1 and t 2 = t 2 at 9 2 = ir/4. 
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